#include "libs.h"
//#include<vector>
//#include<cstring>
//#include<string.h>
#include<iomanip>
#include<SpatialIndex.h>
#include <limits> 
#include<queue>

//#include<cstdlib>
//#include<queue>
//#include<algorithm>
using namespace std; 

extern ISpatialIndex* tree;

const double EPS=1e-9;
const int N=20000;
const int M=51;
const int DIV=100;
const double MAXIMUM=M*DIV;
const int MOD=10007; 

int X[]={-1,0,1,1,1,0,-1,-1};
int Y[]={-1,-1,-1,0,1,1,1,0};

double A[DIMENSION+2][DIMENSION+2]={
  {5.76430151797, 4.38220893224, 0.00433713560189, 7.06315064505},
  {9.90623182304, 2.13645152532, 0.434120821132, 9.05955716438},
  {6.81515602536, 0.824718963702, 8.52522356291, 2.87327829986},
  {1.43620198438, 9.34855160029, 1.36688475402, 3.51637480336}};

double Inverse_A[DIMENSION+2][DIMENSION+2]={
  {-1.1457,   0.7819,   -0.0981,    0.3670},
  {-0.3955,   0.2296,   -0.0507,    0.2444},
  { 0.5087,  -0.3843,    0.1631,   -0.1650},
  { 1.3217,  -0.7803,    0.1114,   -0.4510}};

double T_Inverse_A[DIMENSION+2][DIMENSION+2];
double theta[2*DIMENSION][DIMENSION+2][DIMENSION+2];

node Original_Point[N];
node Encryption_Point[N];
node query_point;

double ServerTime, ProxyTime;
double averagePrecision;
int density_map[M][M];
unsigned int Knn;

struct H_node
{
  int id;
  double dist;
};

struct CMP
{
  bool operator()(const double &A,const double &B)
  {
    return A>B;
  }
};

// implementation of MyVisitor class. 

void MyVisitor::setTarget(Region *t)
{
  target = t;
}

void MyVisitor::visitData(const IData &d) 
{
  //cout << "visiting data ... "  << endl; 
  int id=d.getIdentifier();  
  nodeIndex.push_back(id);
  /*if(validate(id)==true)
    {
      nodeIndex.push_back(id);
      }*/
}

void MyVisitor::visitNode(const INode &n)
{
  //cout << "Visiting node...\n" ;
  visitedNode++;
}

bool MyVisitor::validate(int id) 
{
    int i;
    double ans;
    for(i=0;i<2*DIMENSION;i++)
      {
	ans = Cal_Result(Encryption_Point[id].coordinate,theta[i]);
	if(ans>0)
	  return false;
      }
    return true;
}

void Pre()    //   Get Original data and encrypt original data into encrypted data
{
  FILE *fin=fopen("input1.txt","r");
  FILE *fout=fopen("density_map.txt","w+");
  int i,j,pos[DIMENSION];
  for(i=0;i<DIMENSION+2;i++)
    {
      for(j=0;j<DIMENSION+2;j++)
	{
	  T_Inverse_A[j][i]=Inverse_A[i][j];
	}
    }
  printf("Get the original data from file\n");	
  for(i=0;i<N;i++)   // Get Original Data
    {
      fscanf(fin,"%d%d%d",&(Original_Point[i].operation),&(Original_Point[i].ID),&(Original_Point[i].dimension));
      for(j=0;j<DIMENSION;j++)   // The data in the first d dimension are from raw data
	{          
	  fscanf(fin,"%lf",&(Original_Point[i].coordinate[j]));
	  pos[j]=Convert_Map_Index(Original_Point[i].coordinate[j]);
	}	
      Original_Point[i].coordinate[j++]=1;  // The data in the last two dimension are 1 and random noise
      Original_Point[i].coordinate[j++]=Generate_Random_Noise();
      Original_Point[i].dimension+=2;
      density_map[pos[0]][pos[1]]++;
    }
  printf("start Encryption.\n");
  // Do Encryption for Original Data
  for(i=0;i<N;i++)
    {
      Encryption_Point[i].operation=Original_Point[i].operation;
      Encryption_Point[i].ID=Original_Point[i].ID;
      Encryption_Point[i].dimension=Original_Point[i].dimension;
      Matrix_Multiplication1(Original_Point[i].coordinate,Encryption_Point[i].coordinate);
    }
  for(i=0;i<M;i++)
    {
      for(j=0;j<M;j++)
	{
	  if(j==0)
	    fprintf(fout,"%d",density_map[i][j]);
	  else
	    fprintf(fout," %d",density_map[i][j]);
	}
      fprintf(fout,"\n");
    }
  printf("Encryption is complete\n");
  /*
    Index the Encrypted data
  */
  printf("start Index \n");
  insert();
  printf("Index is complete\n");
  fclose(fin);
  fclose(fout);
}

double Generate_Random_Noise()
{
  int result;
  srand(unsigned(time(NULL)));
  result=rand()%MOD+1;
  return 1.0*result/MOD;
}

void Matrix_Multiplication1(double ma[DIMENSION+2],double ans[DIMENSION+2])
{
  int i,j;	
  for(i=0;i<DIMENSION+2;i++)
    {
      double result=0;
      for(j=0;j<DIMENSION+2;j++)
	result+=A[i][j]*ma[j];
      ans[i]=result;
    }	
}

void Matrix_Multiplication2(double ma1[DIMENSION+2],double
			  ma2[DIMENSION+2],double
			  ans[DIMENSION+2][DIMENSION+2])
{
  int i,j;
  for(i=0;i<DIMENSION+2;i++)
    {
      for(j=0;j<DIMENSION+2;j++)
	{
	  ans[i][j]=ma1[i]*ma2[j];
	}
    }
}

void Matrix_Multiplication3(double ma1[DIMENSION+2][DIMENSION+2],double
			  ma2[DIMENSION+2][DIMENSION+2],double
			  ans[DIMENSION+2][DIMENSION+2])
{
  int i,j,k;
  for(i=0;i<DIMENSION+2;i++)
    {
      for(j=0;j<DIMENSION+2;j++)
	{
	  ans[i][j]=0;
	  for(k=0;k<DIMENSION+2;k++)
	    {
	      ans[i][j]+=ma1[i][k]*ma2[k][j];
	    }
	}
    }
}

void Main()
{
  FILE *fin=fopen("query_point.txt","r");
  int i;
  unsigned int result;
  while(fscanf(fin,"%d",&Knn)!=EOF)
    {      
      fscanf(fin,"%lf%lf",&(query_point.coordinate[0]),&(query_point.coordinate[1]));
      query_point.coordinate[2]=1,query_point.coordinate[3]=Generate_Random_Noise();
      printf("Centered at %lf %lf to find %d nearest neighbor\n",
           query_point.coordinate[0],query_point.coordinate[1],Knn);   
      double high_radius=MAXIMUM,low_radius=0,mid_radius;   
      while(high_radius-low_radius>EPS)
	{
	  mid_radius=low_radius+(high_radius-low_radius)/2;
	  result=Cal(mid_radius);
	  if(result>=Knn)
	    high_radius=mid_radius;
	  else
	    low_radius=mid_radius;
	}
      cout<<"Minimal radius: "<<setprecision(15)<<low_radius<<", number of points: "<<result<<endl;
      printf("Pause, please any number to start\n");          
      scanf("%d",&i);
      /* Based the rough range centered at query point, we use binary
	 search to find the minimal bounding box */
      Binary_Search(low_radius);
    }
  averagePrecision/=100;
  printf("The average precision is %lf\n",averagePrecision);
  printf("The total server time: %lf\n", ServerTime);
  printf("The total proxy time: %lf\n", ProxyTime);
  fclose(fin);
}

void Get_Corner(int len,double plow[DIMENSION+2],double phigh[DIMENSION+2])
{
  int i,j;
  node Original_Corner[4];
  node Encryption_Corner[4];
  Original_Corner[0].coordinate[0]=query_point.coordinate[0]-len;
  Original_Corner[0].coordinate[1]=query_point.coordinate[1]-len;
  Original_Corner[1].coordinate[0]=query_point.coordinate[0]+len;
  Original_Corner[1].coordinate[1]=query_point.coordinate[1]-len;
  Original_Corner[2].coordinate[0]=query_point.coordinate[0]+len;
  Original_Corner[2].coordinate[1]=query_point.coordinate[1]+len;
  Original_Corner[3].coordinate[0]=query_point.coordinate[0]-len;
  Original_Corner[3].coordinate[1]=query_point.coordinate[1]+len;
  for(i=0;i<DIMENSION+2;i++)
    {
      plow[i]=MAXIMUM;
      phigh[i]=-MAXIMUM;
      Original_Corner[i].coordinate[2]=query_point.coordinate[2];
      Original_Corner[i].coordinate[3]=Generate_Random_Noise();
    }
  for(i=0;i<DIMENSION+2;i++)
    {
      Matrix_Multiplication1(Original_Corner[i].coordinate,Encryption_Corner[i].coordinate);
    }
  for(i=0;i<4;i++)
    {
      for(j=0;j<DIMENSION+2;j++)
	{
	  if(Encryption_Corner[i].coordinate[j]<plow[j])
	    plow[j]=Encryption_Corner[i].coordinate[j];
	  if(Encryption_Corner[i].coordinate[j]>phigh[j])
	    phigh[j]=Encryption_Corner[i].coordinate[j];
	}
    }
}

/*
void Proxy_Computing(vector<int> nodeIndex)
{
  FILE *fout=fopen("output.txt","w+");
  priority_queue<H_node, vector<H_node>, CMP> Q;
  int i, size=nodeIndex.size(), id;
  //printf("The number of returned point: %d\n",size);
  double dist;
  for(i=0;i<size;i++)
    {
      id=nodeIndex[i];
      dist=Distance(Original_Point[id].coordinate[0],
      Original_Point[id].coordinate[1]);
      H_node t;
      t.id=id, t.dist=dist;
      Q.push(t);
    }
  fprintf(fout,"\n");
  averagePrecision+=1.0*Knn/nodeIndex.size();
  fprintf(fout,"The precision is %lf\n",1.0*Knn/nodeIndex.size());
  for(i=0;i<Knn;i++)
    {
      fprintf(fout,"%d: %lf\n",Q.top().id,Q.top().dist);
    }
  fprintf(fout,"\n");
}
*/

void Binary_Search(double maxLength)
{
  unsigned int i,j;
  double plow[DIMENSION+2],phigh[DIMENSION+2];
  double low=0,high=ceil(maxLength),mid;
  cout<<"high bound "<<high<<endl;
  int communication_round=0;
  double w[2*DIMENSION][DIMENSION+2]={{0,0,0,-1},{0,0,0,1},{0,0,0,-1},{0,0,0,1}};
  double u[2*DIMENSION][DIMENSION+2]={{1,0,0,0},{1,0,0,0},{0,1,0,0},{0,1,0,0}};
  double uw[2*DIMENSION][DIMENSION+2][DIMENSION+2];
  double temp[2*DIMENSION][DIMENSION+2][DIMENSION+2];
  // Use binary search to find the minimum bounding box that covers K points
  //cout<<"Start Binary Search for Range Query"<<endl;
  time_t t_start=clock();
  while(high-low>EPS)
    {
      mid=low+(high-low)/2;
      Get_Corner(mid,plow,phigh);
      u[0][2]=-(query_point.coordinate[0]-mid),u[1][2]=-(query_point.coordinate[0]+mid);
      u[2][2]=-(query_point.coordinate[1]-mid),u[3][2]=-(query_point.coordinate[1]+mid);
      for(i=0;i<2*DIMENSION;i++)
	Matrix_Multiplication2(u[i],w[i],uw[i]);
      for(i=0;i<2*DIMENSION;i++)
	{
	  Matrix_Multiplication3(T_Inverse_A,uw[i],temp[i]);
	  Matrix_Multiplication3(temp[i],Inverse_A,theta[i]);
	}
      MyVisitor Vis;
      Region r=Region(plow,phigh,DIMENSION+2);
      Vis.setTarget(&r);
      tree->intersectsWithQuery(r,Vis);
      communication_round++;
      cout<<communication_round<<" ends, find "<<Vis.nodeIndex.size()<<" points"<<endl;
      if(Vis.nodeIndex.size()>=Knn)
	high=mid;
      else
	low=mid;
    }
  printf("System Pause, please enter any number to start\n");
  scanf("%d",&i);
  /*
    Based on binary search we can find a minimum_bounding_box that
guarantees to return K nearest neighbors 
  */
  printf("the size of range query %lf\n",low);
  double length=ceil(sqrt(2)*low);
  Get_Corner(length,plow,phigh);
  u[0][2]=-(query_point.coordinate[0]-length),u[1][2]=-(query_point.coordinate[0]+length);
  u[2][2]=-(query_point.coordinate[1]-length),u[3][2]=-(query_point.coordinate[1]+length);
  for(i=0;i<2*DIMENSION;i++)
    Matrix_Multiplication2(u[i],w[i],uw[i]);
  for(i=0;i<2*DIMENSION;i++)
    {
      Matrix_Multiplication3(T_Inverse_A,uw[i],temp[i]);
      Matrix_Multiplication3(temp[i],Inverse_A,theta[i]);
    }
  MyVisitor Vis;
  Region r=Region(plow,phigh,DIMENSION+2);
  Vis.setTarget(&r);
  tree->intersectsWithQuery(r,Vis);
  ServerTime+=(double)(clock()-t_start)/CLOCKS_PER_SEC;
  cout<<"The total number of points:"<<Vis.nodeIndex.size()<<endl;
  t_start=clock();
  //Proxy_Computing(Vis.nodeIndex);
  ProxyTime+=(double)(clock()-t_start)/CLOCKS_PER_SEC;
 
  vector<double> v1;
  vector<double> v2;
  double t_d;
  for(i=0;i<Vis.nodeIndex.size();i++)
    {
      v1.push_back(Distance(Original_Point[Vis.nodeIndex[i]].coordinate[0],Original_Point[Vis.nodeIndex[i]].coordinate[1]));
      //printf(" %lf",Distance(Original_Point[Vis.nodeIndex[i]].coordinate[0],Original_Point[Vis.nodeIndex[i]].coordinate[1]));
    }
  //printf("\n");
	priority_queue<double, vector<double>, CMP> Q;
  for(i=0;i<N;i++)
    {
      t_d=Distance(Original_Point[i].coordinate[0],Original_Point[i].coordinate[1]);
      Q.push(t_d);
    }
  for(i=0;i<Knn;i++)
    {
      v2.push_back(Q.top());
      //printf(" %lf",Q.top());
      Q.pop();
    }
  //  printf("\n");
  for(i=0;i<v2.size();i++)
    {
      for(j=0;j<v1.size();j++)
	{
	  if(v1[j]==v2[i])
	    break;
	}
      if(j==v1.size())
	break;
    }
  if(i<v2.size())
    {
      printf("Incorrect\n");
      printf("System Pause, Please enter any number to start\n");
      scanf("%d",&i);
    }
  else
    printf("Correct\n");
}

int Cal(double radius)
{
  queue<cell> Q;
  int i,j;
  bool vis[M][M];
  int pos_x,pos_y;
  for(i=0;i<M;i++)
    {
      for(j=0;j<M;j++)
	{
	  vis[i][j]=false;
	}
    }
  struct cell now;
  struct cell next;
  int sum_point=0;
  pos_x=Convert_Map_Index(query_point.coordinate[0]);
  pos_y=Convert_Map_Index(query_point.coordinate[1]);
  if(Check_In_Circle(pos_x,pos_y,radius)==true)
    {
      now.x=pos_x;
      now.y=pos_y;
      Q.push(now);
      vis[now.x][now.y]=true;
    }
  while(!Q.empty())
    {
      now=Q.front();
      sum_point+=density_map[now.x][now.y];
      Q.pop();
      for(i=0;i<8;i++)
	{
	  pos_x=now.x+X[i];
	  pos_y=now.y+Y[i];
	  if(pos_x>=0&&pos_x<M&&pos_y>=0&&pos_y<M)
	    {
	      if(vis[pos_x][pos_y]==false&&Check_In_Circle(pos_x,pos_y,radius)==true)
		{
		  vis[pos_x][pos_y]=true;
		  next.x=pos_x;
		  next.y=pos_y;
		  Q.push(next);
		}
	     }
	}
    }
  //cout<<"The number of total point is "<<sum_point<<" in the radius of "<<setprecision(15)<<radius<<endl;
  return sum_point;
}

bool Check_In_Circle(int now_x,int now_y,double radius)
{
  int i;
  int corner_x[4],corner_y[4];
  corner_x[0]=DIV*now_x,corner_x[1]=DIV*(now_x+1),corner_x[2]=corner_x[1],corner_x[3]=corner_x[0];
  corner_y[0]=DIV*now_y,corner_y[1]=corner_y[0],corner_y[2]=DIV*(now_y+1),corner_y[3]=corner_y[2];
  for(i=0;i<4;i++)
    {
      if(Distance(corner_x[i],corner_y[i])+EPS>radius)
	return false;
    }
  return true;
}

double Distance(double x,double y)
{
  return sqrt((x-query_point.coordinate[0])*(x-query_point.coordinate[0])+(y-query_point.coordinate[1])*(y-query_point.coordinate[1]));
}

int Convert_Map_Index(double num)
{
  int ans=num/DIV;
  return ans;
}

void Main1()
{
  int i,j;
  FILE *fin1=fopen("test_query.txt","r");
  FILE *fin2=fopen("input1.txt","r");
  double plow[2],phigh[2],t;
  while(fscanf(fin1,"%lf%lf%lf%lf",plow,plow+1,phigh,phigh+1)!=EOF)
    {
      printf("%lf %lf\n",plow[0],phigh[0]);
      printf("%lf %lf\n",plow[1],phigh[1]);
      MyVisitor Vis;
      Region r=Region(plow,phigh,2);
      Vis.setTarget(&r);
      tree->intersectsWithQuery(r,Vis);

      vector<int> Point;
      for(i=0;i<N;i++)
	{
	  for(j=0;j<2;j++)
	    {
	      fscanf(fin2,"%lf",&t);
	      if(t<plow[j]||t>phigh[j])
		break;
	    }
	  if(j==2)
	    Point.push_back(i);
	}
      printf("Point:%d, nodeIndex:%d\n",Point.size(),Vis.nodeIndex.size());
      if(Point==Vis.nodeIndex)
	printf("yes\n");
      else
	printf("no\n");
    }
}
